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An Accuracy Evaluation of Unstructured 
Node-Centred Finite Volume Methods 

Magnus Svard* Jing GongUnd Jan Nordstrom,* 


Abstract 

Node-centred edge-based finite volume approximations are very 
common in computational fluid dynamics since they are assumed to 
run on structured, unstructured and even on mixed grids. 

We analyse the accuracy properties of both first and second deriva- 
tive approximations and conclude that these schemes can not be used 
on arbitrary grids as is often assumed. For the Euler equations first- 
order accuracy can be obtained if care is taken when constructing 
the grid. For the Navier-Stokes equations, the grid restrictions are 
so severe that these finite volume schemes have little advantage over 
structured finite difference schemes. Our theoretical results are veri- 
fied through extensive computations. 


1 Introduction 

Finite volume approximations are widely used in computational fluid dy- 
namics (CFD). There are several different variations of finite volume ap- 
proximations and one of the more popular is the node-centred edge-based 
approximation. (See [1-13].) One reason for its popularity is the simple 
data structures associated with this scheme which make aerodynamic com- 
putations very efficient. Another is that it is assumed to run on grids made 
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Figure 1: A typical CFD finite volume grid for a wing. 


up of any type of elements, such as quadrilaterals or triangles in two space 
dimensions or tetrahedrons, prisms etc in three dimensions. Not only is the 
scheme assumed to yield at least first-order accuracy on all these element 
types but it is also assumed that a grid may mix element types. Such a grid 
is seen in Fig. 1. (See [14] and also [13, 15, 16] for more hybrid grids.) This 
property was called grid transparency in [1] and it is a crucial property since 
it greatly simplifies the task of grid generation. 

In a sequence of previous articles [8, 9, 17], the present authors have 
investigated the accuracy and stability of the edge-based finite volume ap- 
proximation for both hyperbolic and parabolic problems. In this paper we 
focus on accuracy questions for Cauchy problems. We review and extend the 
analysis of the previous papers to give a coherent view of the accuracy of the 
edge-based finite volume approximation. 

Often, it is assumed that by only including the Laplacian part of the 
viscous terms, a good approximation for the Navier-Stokes equations is ob- 
tained and that makes it possible to construct a compact discretisation for 
the viscous terms. 

We will consider the scalar advection-diffusion equation as a model for 
the Navier-Stokes equations, 

u t + au x + bu y = e(u xx + v yy ). (1) 

This equation contains the most important terms appearing in the Navier- 
Stokes equations and if (1) is accurately approximated, the same is true for 
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the Navier-Stokes equations. 

The contents of this report are divided as follows; in Section 2 we will 
present the edge-based finite volume approximations; then follow Section 3 
and 4 containing accuracy studies for the Laplacian approximation and the 
advection equation respectively. Finally, conclusions are drawn in Section 5. 


2 The Finite Volume Approximation 


We begin by stressing that neither of the schemes derived below are due 
to the present authors but rather standard schemes used in CFD (c.f. [1]- 
[13], [17], [18]). 

Following the derivation in [8], we begin with one advective term and 
consider equation (1) with b = e = 0. The finite volume approximation is 
derived from the weak form of the equation, 


IL 


u t dxdy 



(2) 


Before we proceed to approximate (2) we introduce some notation. The 
discrete solution will be defined at a grid vertex (c.f Fig. 2). (That defines 
a node-centred scheme.) Let r* denote a grid point and let V) be an n- 
sided polygon with sides dsi n . V) is defined as the volume inside the dual 
grid around r,. The dual grid is in turn defined as the straight lines drawn 
between the centres of mass of the cells with r t as a vertex and the midpoints 
of the edges from r,, see Fig. 2. Further, dsi n is defined as the sum of the 
length of the “centre of mass-midpoint-centre of mass” lines passing over one 
edge (see Fig. 2). We will denote the measure of the volume V t , although it 
is a slight abuse of notation since the volume itself is also denoted V). 

Further, let (wjv)m denote the outward pointing derivative normal to dsi n . 
Let ri n = \r\ — r n \. Finally, let iV, denote the set of indices of points being 
neighbours to r,. 
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Figure 2: A generic 2D grid. Solid lines are the grid lines and dashed lines 
corresponds to the dual grid. 


Now we return to equation (2). The integration is carried out over a dual 
volume of size V % for point r* (Fig. 2). 



udy = 0 


( 3 ) 


Along each dsi n we approximate the solution by (u n + iq)/2, where u n and 
Ui are the discrete solution values at r n and r,. We obtain, 


~\ r ( \ i Ui + U n A n 

Vi(ui)t + a 2 _^ — 2 — Ayn = 

nENi 


( 4 ) 


where the sum ranges over all neighbours to and A y n is the difference in 
y along dsi n . A similar approximation is obtained for the y-derivative. 

On a rectangular grid with a smooth mapping to a Cartesian grid, the 
first derivative approximation in (4) can be proven second-order accurate 
using Taylor expansions. It is also possible to prove first-order accuracy on 
an unstructured triangular grid (see Appendix I). This is an upper bound of 
the error and on highly regular grids the error may even be second-order. 

As mentioned in the introduction, a common approximative model for 
the viscous terms used in computations with the Navier-Stokes equations is 
to only include the Laplacian. Following the derivation in [9], assume that 
a = b = 0 and e = 1 in equation (1) and integrate over the domain Vj. We 
obtain, 
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where Gauss’ theorem is used. N denotes the outward pointing unit normal 
vector such that J^- = = Vu-N. A straightforward approximation of (5) 

would be, 

^ ^ (uN^indSin- (6) 

ndNi 

The normal derivative to dsi n (('Ujv)m), is approximated by the central dif- 
ference along that edge and we obtain for an interior point r*, 

Vi(u t )i = ^2 — — — ds in- (7) 

_ Ar I'm 

neNi 

Using Taylor expansions this approximation can be shown second-order ac- 
curate on Cartesian grids. 


3 Accuracy of Laplacian on unstructured grids 

In [9] the approximation (7) of the Laplacian was studied and a stable imple- 
mentation of the boundary conditions derived and several computations were 
made for the wave equation u tt = A u in two space dimensions. A grid con- 
vergence study was made and the results are shown in Fig. 3 for a Cartesian 
mesh. Second-order accuracy is obtained in this case. 



Figure 3: Grid convergence in l 2 , for the wave equation on a Cartesian mesh. 
The dashed reference line represents second-order of accuracy. 

Next, we tested the scheme on a triangulated Cartesian grid (see Fig. 4). 
As this kind of grid is refined, the convergence results are poor as seen in Fig. 
4. The reason for the bad grid convergence will be analysed in the following 
subsection. 
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Figure 4: Left: A triangulated Cartesian mesh. Right: ^-convergence on a 
triangulated Cartesian mesh. 

Remark A triangulated Cartesian grid is not truly an unstructured grid. 
But it is a good model of an unstructured grid and if the scheme is incon- 
sistent on such a grid it will not be consistent on an arbitrary unstructured 
grid. 

3.1 Consistency Analysis 

We will analyse consistency for the centre point of the grid in Fig. 5. The 
scheme (7) leads to, 

Vi(ui)t = ^2 ~ — ~ ds ^ ( 8 ) 

i 2 Til 

where V\ = K 2 is the area inside the dashed line. By Taylor expansions we 
obtain for the edge r 2 1- 

u 2 ~ ui _ «i + hu x + huy + \{h 2 u xx + 2 h 2 u xy + h 2 u yy ) - + 0(h 3 ) 

V2h ~ y/2 h 

where the derivatives are taken at point 1. A similar derivation for edge r 5i 
leads to, 

U2 , U5 xi\ , h u xx T 2 h u xy -[- hi u yy -f- &(^h ) 

— (IS 19 1 7 = (IS 15 = . 

V 2 h x/2 h 3 

In the same manner we obtain for points edges r 3 i,r4i,r 6 i and rn- 

— ^ — ds 13 H — ds i6 — u xx + 0 {h )), 

u 7~u± U4 — Ui _ VS fh2 . //-)/l3\\ 

— - — ds 17 H — dsi4 — u yy + Oyh )). 
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Summing up, using (8) and V\ = K 2 yield, 

/ \ Id - y/5 2 . 

(«i )t = — g — A m + + 0(h), 

i.e. an 0(1) error. This is stated in theorem 3.1. 

Theorem 3.1 On a general grid the approximation (1) of (1), with a = b = 
0 and e = l, is inconsistent. 
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Figure 5: A triangular grid where the dual grid is the dashed line, b = h\/5/3 
and a = h\/ 2/3. 

The theorem does not imply inconsistency on all grids. It is easily seen in 
the above analysis that if the scheme consists of equilateral polygons, consis- 
tency is recovered. This is verified with computations on a grid with equilat- 
eral triangles. Again, the wave equation is considered and a grid refinement 
is performed. In this case the domain as well as the cells are equilateral 
triangles and the grid convergence corroborates second-order accuracy. (See 
Fig. 6.) 
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Figure 6: Left: A mesh with equilateral triangles. Right: Convergence for the 
wave equation with mesh consisting of equilateral triangles. The reference 
line represents second-order of accuracy. 

It has been shown that the scheme works satisfactory on Cartesian as well 
as equilateral triangular grids which are both equilateral polygons. In fact, 
the scheme will work on grids consisting of any type of equilateral polygons. 
In the case of equilateral triangles the domain was also chosen to be a triangle. 
It is not possible to alter the shape of the cell triangles near the boundary 
since that would introduce errors of order 1. 

The triangulated Cartesian grid used in Fig. 5 resulted in an inconsis- 
tent scheme. In the equilateral triangular case the scheme coincides with a 
mass lumped finite element scheme. To show the differences between finite 
volume and finite element schemes we derive the mass lumped finite element 
approximation of Laplace equation, with linear basis functions applied to the 
grid in Fig. 5 and obtain, 


h 2 (ui) t = -v A ~v 3 + 4ui - u 7 - u 6 . 


The left hand side, h 2 (ui) t , is the result of the mass lumping and is precisely 
equal to obtained with the finite volume technique. However, the 

right hand side coincides with the ordinary finite difference scheme on the 
corresponding Cartesian grid. Compared to the finite volume case, different 
weights on all the neighbouring points are obtained and especially the cross 
derivative contribution is cancelled since u 5 and u- 2 do not enter the scheme. 

Remark We use the wave equation in our numerical computations, since 
it is a less forgiving equation. However, also for the heat equation on a 
triangulated Cartesian grid, the convergence levels out (see Fig. 7). 
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Figure 7: Convergence for the heat equation on a triangulated Cartesian 
grid. 

3.2 A Different Finite Volume Approximation 

To further explore the finite volume discretisations of viscous terms we will 
briefly discuss a discretisation proposed in [7] and also used in [6]. This 
discretisation allows for the computation of all viscous terms in the Navier- 
Stokes equations but it is not compact as the Laplace approximation previ- 
ously discussed. The scheme is basically an application of the edge-based first 
derivative approximation (4) twice. In [7] it is observed that this scheme does 
not take into account the closest neighbours possibly yielding poor smoothing 
properties. (Compare the second derivative approximation that is obtained 
by applying a centred finite difference scheme twice.) Accordingly, they pro- 
pose an augmented form that also uses the neighbours. 

We will analyse consistency of the first derivative applied twice without 
this augmentation to see whether this can be used as an approximation of 
the Laplacian on more general grids. It certainly has the disadvantage that 
it is wider and thus not as computationally efficient. Recalling that N t is the 
neighbours of the point indexed i and denoting by (u x )i the approximation 
of the x-derivative at the same point we have that, 

K)i = y; — 2 — A% ’’ 

’ j£Ni 

where A i/j is the difference (with sign) in the y direction going anti-clockwise 
along dsij. If approximation (9) is carried out for point 1 in Fig. 5, we would 
obtain with V\ = li 2 , 

(«*) i = y ^ 2 ~ Ayj = Ux ^ Vl ^ + ( 10 ) 

j£N i 
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Further, if the grid in Fig. 5 is extended in all directions we obtain second 
order approximations similar to equation (10) of ( u x \ at all the points 2 — 7. 
Next, applying (9) for point 1 again, yields a first-order approximation of 

xx)i • 

l \ _ 1 ( Ux )j + ( U x)l + 0(^ 2 ) A _ / \ x 

i'^xx) i — Ty / J ~ ^Vj — yi)xx H - (11) 

V i z J a 

jeN i 

However, if the grid is not extended precisely as in Fig. 5, for example 
if the location of the points 2 and 5 are slightly disturbed or perhaps the 
diagonal between 1 and 5 is replaced by a diagonal between 4 and 6, a term 
of the form hu xy will enter the Taylor expansion of (10). As an example, 
consider point 1 with the proposed diagonal altered. (Compare also to Fig. 
4 where the orientaion of the diagonals are somewhat arbitrary.). Then, 

(u x ) i = u x (x i, yfi) + const ■ hu xy + 0(h 2 ). (12) 

Hence, we have shown that a first-order error may occur in the first derivative 
approximation. Then, suppose that (u x ) 5 has an error of the form (12). Then 
( u x )i , i G {2, 3, 4, 6, 7} in Fig. 5 are approximated to second-order and (u x ) 5 
only to first-order. Applying (9) again to compute (u xx ) 1 yields, 

j N 1 ( Ux )j + ( U x)l + 0(h) A _ / \ y nn s / 1q \ 

y^xx)l -y / J 2 ^Vj U\X \ ? Vijxx ^(l)* (13) 

1 jeN 1 

Note that since (u x ) 5 = u x (x 5 , y. 5 ) + 0(h) is the only first order term there 
is no other term that can cancel that error. Multiplying by ~ yields 
an error of 0(1). This example shows that unless the grid has a high degree 
of regularity the scheme becomes inconsistent. The above considerations can 
be summarised in the following theorem. 

Theorem 3.2 On a general grid, the application of the first derivative ap- 
proximation (9) twice to approximate a second derivative results in an incon- 
sistent scheme. 

Remark Note that the above theorem does not imply that the scheme is 
inconsistent on all grids. As shown above, with certain restrictions imposed 
on the grid, the scheme is consistent. 

As mentioned above, in [7] they propose an augmented form of the above 
approximation. However, the augmentation assumes implicitly that the 
above approximation is consistent. Thus, the same restriction of possible 
grids applies to their scheme as well. 
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4 The advection equation on mixed grids 

We have already mentioned that the first derivative approximations do result 
in at least first-order accuracy on unstructured triangular grids (see [8] and 
Appendix I). However, as the previous section reveals, the finite volume 
approximations are not consistent on all grids, though they are assumed to be 
in real life computations. Hence, we go on studying a common technique used 
in CFD when constructing grids. That is to use structured Cartesian grids in 
boundary layers which is changed to triangular grids outside boundary layers 
which more easily adapts to complex geometries. Throughout, this section 
we will assume that e = 0 in equation (1). 

4.1 Consistency analysis of interface 

In this case we begin by analysing the first derivative approximation at an 
interface between unstructured triangles and quadrilaterals. In particular, we 
consider the finite volume approximation at an interface between arbitrary 
unstructured triangles triangular and a Cartesian grid. (See Fig. 8) 



Figure 8: An interface between triangles and squares. The dashed line marks 
the dual grid. 

We assume that the side of the squares have length h. This interface 
is smooth and we can derive the expected order of accuracy at the point 
c. The standard way of deriving order of accuracy is to make use of Taylor 
expansions. However, that becomes very complicated when unstructured 
triangles are considered. Instead, we will consider the two approximations 
obtained if the grid is split along the interface. Then we formally need to 
specify u c on both sides using boundary conditions, but in this case u c is 
given by the other side. First consider the finite volume approximation on 
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the control volume related to 1 — 2 — c — 5. Denote the corresponding control 
volume Vf? . (Note that the boundary of V{r partly coincides with 2 — c — 5) 
The approximation (4) becomes, 

v c ( «c)t = Ul + Uc h + u c (-h ) (14) 

where c is a boundary point and the flux through the boundary is approx- 
imated by u c . This is a consistent treatment of the boundary (see [8]) and 
the scheme is consistent. 

The finite volume approximation (4) on the triangles becomes, 

V*(u c ) t = ^2 Ul+ 2 Uc ^Vi + u c h, (15) 

%= 2 

where V* is the corresponding control volume related to c — 2 — 3 — 4 — 5. 
(Again, it partly coincides with 2 — c— 5 where it precisely matches V c c . ) This 
is also a consistent approximation according to [8]. Note that, + Vf? = V c . 
Hence, the sum of the two approximations is, 




Ui + u c 


Ui + Ur 


A</j. 


(16) 


i = 2 


i = 1 


Note that, (16) is the same finite volume approximation at u c as if the scheme 
had been applied directly ignoring the interface. Since it is constructed from 
two first-order approximations the sum is also at least first-order. 

Since the finite volume approximation on a smooth quadrilaleral grid is 
also second-order accurate we conclude that the above reasoning applies to 
an interface between any smooth equilateral grid and triangles. 

In order to derive the global order of accuracy we will use some theory 
from [8]. Consider the equation u t = u x + F(x ) with a finite volume dis- 
cretisation v t = Dv + F where v denotes the vector with components Vi. D 
is the resulting matrix obtained by the finite volume discretisation, and F 
is a forcing function. According to [8] there is an energy estimate such that 
||u|| < constant^ ||/|| + ||g|| + ||T||) where f,g denotes the initial data and 
boundary data. || • || is the discrete / 2 -norm (in two space dimensions) defined 
by, 




(17) 


Then the error equation is, 


e t = De + T, 


(18) 
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where e — u — v and T is the truncation error. Note the e has zero initial 
and boundary data. 

To derive the global order of accuracy we consider a hybrid grid with N 2 
points in two space dimensions and a typical grid size 0(h) = 0(1/N). The 
number of interface points is 0(N) and V % ~ K 2 . Then from (18) we have 
||e|| < constant\\T\\, i.e. the size of ||e|| is proportional to the size of ||T||. If 
we assume that the triangular grid is highly symmetric such that the local 
order of accuracy is 2 and the order of accuracy at the interface is 1 (as shown 
above) we obtain, 

||T|| 2 = ^T 2 Vi ~ (O(N)O(h 2 ) + 0{N 2 )0{h i ))0{h 2 ) = 0(h 3 ). (19) 

Hence ||e|| Ill’ll G(h 15 ), i.e. the order of accuracy is 1.5. (If the accuracy 
on the triangular part is only first order the global order order of accuracy 
will be 1.) 

Next, consider the interface in Fig. 9. 



Figure 9: An interface between triangles and squares. The dashed line marks 
the dual grid. 

If we assume that a grid refinement would keep the corner we obtain from 
Taylor expansions with the sides of the squares have length h. Then V c = | K 2 
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and, 


h 2 

U\ — u c + hu x T ~^~u xx T 0(lf ) ■ 
ti 2 

U'2 — u c hiiy T ^ ^yy T 0(h ), 
h 2 

1/3 U c hUx hlly T ~^~{Uyy T U XX “I" 2 U X y) T O (ll ), 

h 2 

M4 — u c hu x -t- ~u xx -I - 0(h ), 
h 2 

u 5 — u c + htiy + ^ ^yy T ^ ^ ) - 


The finite volume approximation (4) becomes, 


u c + u 1 u c + u 2 h u c + u 3 (-h) u c + u 4 (-5h) 

^ + - 2-6 + — 3 " + — 6 “ = 

K 2 

V c u x (x c , y c ) + — u y (x c , y c ) + (9(fi 3 ). (20) 


Shifting the diagonal in Fig. 9 will not improve things. In either case, we 
have an 0(1) error at c. 

If we assume that we have an interface with a number of corners like the 
one above. With two space dimensions and a total number of points 0(N 2 ) 
there are 0(N) corner points with an 0(1) error. Again, we consider the 
error equation and use the energy estimate derived in [8] to conclude that 
||e|| ~ ||X||. To determine the size of the truncation error, assume that the 
order of accuracy away from corners is 1, i.e T t ~ O(h) = 0(1 /N) at those 
points. 


||T|| 2 = ^T 2 h 2 ~ 0(h 2 )(0(N)0(l) 2 + 0(N 2 )0(1/N) 2 ) = 0(h). (21) 

i 


Then ||T|| ~ 0(h 1 ! 2 ) ~ e, i.e the order of accuracy is 0.5. 

Theorem 4.1 Consider a hybrid mesh composed of unstructured triangles 
and smooth quadrilaterals with a non-smooth interface. Then the approxi- 
mation (4) of (1) with b = e = 0 is inconsistent at the interface and the 
global order of accuracy reduces to 0.5. 

Note that throughout this article we consider a model problem with two 
space dimensions and we assume implicitly in Theorem 4.1 that an interface 
is one-dimensional. In general we consider an interface to be a d— 1 surface in 
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Q 

Ntot 

Qi 

N b 

q t ot 

interface type 

2 

N 2 

1 

N 

1.5 

smooth, 2D 

1 

N 2 

1 

N 

1 

smooth, 2D 

2 

N 2 

0 

N 

0.5 

non-smooth, line, 2D 

2 

N 2 

0 

1 

1 

non-smooth, point, 2D 

2 

N 3 

1 

N 2 

1.5 

smooth, 3D 

1 

N 3 

1 

N 2 

1 

smooth, 3D 

2 

N 3 

0 

N 2 

0.5 

non-smooth, plane, 3D 

2 

N 3 

0 

N 

1 

non-smooth, line, 3D 

2 

N 3 

0 

1 

1.5 

non-smooth, point, 3D 


Table 1: Order of accuracy for hybrid meshes, q is the lowest order of 
accuracy of the scheme away from the interface. N tot is the total number 
of points. (We assume N points in one space dimension.) q t is the order of 
accuracy at interface points. (0 meaning order 1 error.) N b is the number of 
boundary points. q to t is the resulting overall order of accuracy. 


a problem with d space dimensions. One can anticipate even more restrictive 
constraints in higher-dimensional problems. If there is a disconituity at one 
or many points or even along a subsapce of the interface order 1 errors will 
be introduced. With the the same technique as used above to estimate the 
/ 2 -norm of the truncation error we get the results displayed in Table 1 in two 
and three space dimensions. 

4.2 Computations 

Consider, equation (1) with a = 1, b = 1 and e = 0 on a two-dimensional rect- 
angular domain discretised by (4). The initial data is u(x , 0) = asin{^{ax + 
by)) and the exact solution, u(x, t) = sin(4:ii(a(x — t) + b(y — t))). The exact 
solution is also used as boundary data. 

In Fig. 10 the convergence rate for a triangular grid is displayed. The 
convergence rate is close to 2 which is due to the fact that the grid is very 
symmetric. 
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Figure 10: Left: A triangular mesh with 727 grid points. Right: Convergence 
rate for the advection equation on triangular meshes. The dashed reference 
line represents second-order of accuracy. 

We also consider a quadrilateral grid (see Figure 11) where there is a 
smooth mapping to a Cartesian grid. In this case we also obtain second- 
order accuracy. See Table 2 



Figure 11: A smooth quadrilateral grid. 

To illustrate that the smoothness is essential for quadrilateral grids we 
show the convergence properties for quadrilateral grids (without stretching) 
with a 20% perturbation in the x- and y-direction for each node. This yields 
a sequence of non-smooth quadrilateral grids. (See Figure 12). 
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Vn 

^-error 

^-convergence 

l 0 0 — error 

41 

0.0325 

- 

0.2497 

81 

0.0080 

2.06 

0.1054 

121 

0.0035 

2.06 

0.0415 

161 

0.0020 

1.98 

0.02195 

201 

0.0013 

1.92 

0.0152 

241 

9.05e-4 

2.00 

0.0107 

281 

6.76e-4 

1.90 

0.0075 


Table 2: Errors and convergence on smooth quadrilateral grids. 


Vn 

^-error 

^-convergence 

loo — error 

41 

0.0424 

- 

0.1187 

81 

0.0225 

0.93 

0.1243 

121 

0.0192 

0.40 

0.1393 

161 

0.0179 

0.25 

0.1534 

201 

0.0168 

0.29 

0.1551 

241 

0.0161 

0.23 

0.1382 

281 

0.0160 

0.04 

0.1855 


Table 3: Errors and convergence on non-smooth quadrilateral grids. 



Figure 12: A non-smooth quadrilateral grid. 

In Table 3 it is seen that there is no convergence in l ^ and in l 2 the 
convergence drops towards 0. The results in shows that there are 0(1) 
errors present at at least one point. Since the / 2 -error does not go to zero we 
conclude that the number of inconsistent points are proportional to the total 
number of points. This is what we would expect since there are perturbations 
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Vn 

/ 2 -error 

^-convergence 

loc — error 

41 

0.0352 

- 

0.2603 

81 

0.0087 

2.05 

0.1035 

121 

0.0038 

2.06 

0.0412 

161 

0.0022 

1.91 

0.0353 

201 

0.0014 

2.04 

0.0291 

241 

0.0010 

1.85 

0.0246 

281 

7.5e-4 

1.87 

0.0218 


Table 4: Hybrid quadrilateral and triangular grid with smooth interface. 


to the nodes everywhere. 

Next, we turn to mixed grids. We begin our study with the same smooth 
quadrilateral grid used above. The total number of points is (N + l) 2 and 
we place an interface at x(N/2 + 1). This yields a smooth interface. On 
one side we triangulate the grid and on the other we keep the quadrilaterals 
(See Figure 13). According to the previous theory this should have order of 
accuacy 1.5. The errors and convergence are displayed in Table 4 and we see 
that the errors in l 2 and /<*, are convergent. 



Figure 13: A hybrid quadrilateral and triangular grid with a smooth interface. 

However, as previously mentioned, it is common that the interface is 
non-smooth. (See Figure 1.) An example of such an interface is displayed in 
Figure 14 and the convergence data given in Table 5. 
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Vn 

Z 2 -error 

^-convergence 

l 0 0 — error 

41 

0.0338 

- 

0.1329 

81 

0.0144 

1.23 

0.1015 

121 

0.0108 

0.72 

0.1003 

161 

0.0091 

0.60 

0.0999 

201 

0.0081 

0.52 

0.1004 

241 

0.0074 

0.50 

0.1006 

281 

0.0068 

0.55 

0.1004 

321 

0.0064 

0.46 

0.1004 


Table 5: Hybrid quadrilateral and triangular grid with non-smooth interface. 



Figure 14: A hybrid quadrilateral and triangular grid with a non-smooth 
interface. 

Clearly, there are points that are non-convergent according to Table 5. As 
shown above the ^-convergence should approach 0.5 which is corroborated 
in Table 5. 


5 Conclusions 

Edge-based node-centred finite volume approximations are widely used in 
computational aerodynamics for one important reason. They are assumed 
accurate on unstructured grids. 

In Theorem 3.1 and 3.2, we prove that discretisations of second deriva- 
tives with two different commonly used approximations are inconsistent on 
unstructured grids. However, the accuracy is recovered on grids with a high 
level of regularity. For the Laplacian approximation equilateral polygons 
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are required and for the first derivative approximation applied twice slightly 
more general grids can be allowed. 

The CFD community (see for example [18]) have recognised that in 
boundary layers, Cartesian grids should be used in order to obtain better 
accuracy (see Fig 1). Our analysis confirm that observation since the viscous 
terms are dominant in the boundary layer and those are only approximated 
correctly on regular grids, such as rectangles and equilateral polygons. 

Structured grids in boundary layers and unstructured to capture geome- 
tries have led us to study mixed grids for the first derivative approximation. 
That approximation is shown to be consistent on triangular grids and on 
smooth quadrilateral grids. For mixed grids with smooth interfaces between 
the quadrilaterals and triangles the accuracy is not degraded. 

However, if the interface is non-smooth the order of accuracy drops to 
0.5 due to inconsistencies in the discretisation. To summarise, care has to 
be taken when constructing hybrid grids. Extensive numerical experiments 
corroborate our theoretical results. 
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APPENDIX 


I First derivative approximation 

We will show that the first derivative approximation is first order accurate 
on a general triangular grid. This derivation is due to [19]. 

Consider the grid in Fig. 15. 



Figure 15: A general triangular grid. 


The approximation of u x at the centre point c is, 

/ \ 1 x u c T Ui . . 

Me = yYl — 2 — Ayi: ( 22 ) 

c i 

where V is the volume the measure of the volume of the dual grid. If this 
approximation is first order accurate it should exactly differentiate a constant 
and a linear function. 

For a constant function u = C we have u x = 0 and, 

Me = y ^2 Ayi = ( 23 ) 


Hence, the approximation is correct for a constant. Next, we turn to a linear 
function u = ax + by where a and b are constants and u x = a. Denote by 
?i = (xi,yi) the position of the centre of mass for the triangle 1,2, c and f 2 
for the triangle 2, 3, c etc. Then the approximation for the linear function is, 


1 s~^ax c + by c + axi + byi _ . 

Me =77 2^ 5 W - Vi- 1) = 


^ max . 7 

1 axi + byi . „ v 

— 5 — (Vi-Vi- 1), 
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where yo is interpreted as y ma x in a cyclic manner around c. (In Fig 15 
max = 5.) Hence, it is always assumed assumed that we compute modulo 
the number of neighbours. 

The y-coordinate of the centre of mass is obtained by, 


Vi + Vi+i + He 

Vi = ^ 


(25) 


Then, 


Vi ~ Vi - 1 


Vi + Vi + 1 +y c Vi + Vi+1 + Uc Vi+1 - Vi-1 


(26) 


Using (26) in (24) results in, 

, x _ 1 ^ ax i y i+ 1 - j/i -1 byj y i+1 - y^ 1 

[U x ) c ~ y 2^ 2 3 ' 2 3 

c i 


Note that, 


Hence, 


max 

- Vi-i) = 0 . 

i= 1 


( ) c 


a 

V c 


E Vi+l ~ Vi -1 
Xi 6 


(27) 


Next, we will compute V c . The triangle 1, 2, c has the area, 

\\( r i~r c ) x (r 2 — r c ) | 

where r* = (x*, yi) is the point vector at point i. With the same notation as 
for fi we denote the area of traingle i, i + 1, c by A t . Hence, 

A = ^((xi - x c )(y i+ 1 - y c ) - (y t - y c )(x i+1 - x c )). 

The area of the dual volume is, 

V i = ^2 = S ~ X c)(Vi + 1 “ Vc) ~ (Vi ~ Vc)(Xi+l ~ x c )) = 

i i 

^y^ y Xiy i+1 - Xiy c - x c y i+ i + x c y c - ( yix i+1 - y c x i+1 - y { x c + y c x c ) = 

i 

^ - ViXi+i = ^ ^2xi(y i+ i - yi-i) (28) 

« i 
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Note that we use modulo calculations extensively which also justifies the shift 
of the indices. With (28) in (27) we have, 

(^x)c — ev- 
idence, we have shown first order accuracy on an arbitrary triangular grid. 
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